Mon. Not. R. Astron. Soc. 000, [TH22] (2011) Printed 30 September 2011 (MN WT^ style file v2.2) 

Modelling neutral hydrogen in galaxies using cosmological 
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ABSTRACT 

The characterisation of the atomic and molecular hydrogen content of high-redshift 
galaxies is a major observational challenge that will be addressed over the coming 
years with a new generation of radio telescopes. We investigate this important issue by 
considering the states of hydrogen across a range of structures within high-resolution 
cosmological hydrodynamical simulations. Additionally, our simulations allow us to in- 
vestigate the sensitivity of our results to numerical resolution and to sub-grid baryonic 
physics (especially feedback from supernovae and active galactic nuclei) . We find that 
the most significant uncertainty in modelling the neutral hydrogen distribution arises 
from our need to model a self-shielding correction in moderate density regions. Future 
simulations incorporating radiative transfer schemes will be vital to improve on our 
empirical self-shielding threshold. Irrespective of the exact nature of the threshold we 
find that while the atomic hydrogen mass function evolves only mildly from rcdshift 
two to zero, the molecular hydrogen mass function increases with increasing redshift, 
especially at the high-mass end. Interestingly, the weak evolution of the neutral hy- 
drogen mass function is insensitive to the feedback scheme utilised, but the opposite 
is true for the molecular gas, which is more closely associated with the star formation 
in the simulations. 

Key words: galaxies: evolution: general - galaxies: luminosity function, mass func- 
tion - dark matter - radio lines: galaxies - methods: N-body simulations 



1 INTRODUCTION 



Approximately 75 per cent of the baryonic material in the 
Universe is hydrogen, so understanding and being able to 
reproduce the distribution of its various phases is a crucial 
step in constructing a complete theory of galaxy formation. 
After re-ionisation completed at 2; ~ 6, the majority of the 
hydrogen is expected to reside in a diffuse, ionised state out- 
side of galaxies. However, it is also believed that significant 
reservoirs of neutral (atomic and molecular) hydrogen must 
exist, particularly at high-redshift, as this is the basic fuel 
for star formation. 

At higher densities than the diffuse, ionised intergalac- 
tic medium surrounding the galaxies, large-scale gas clouds 
are visible in absorption as the Lyman-a forest. The latter 
can be probed by observations of intervening low column 
density HI structures in the spectra of quas ars found, for ex- 
ampl e, in the Sloan Digital Sky Survey (e.g. lProchaska et al.l 
120101 ). As we probe ever closer to the galaxy disks the gas 



density increases, along with the HI surface density, allow- 
ing us to directly probe the emission from neutral hydro- 
gen. While a challenging observation, it greatly expands 
our understanding of galaxy formation on scales at which 
the galaxies themselves reside. For these reasons, great ef- 
forts have been made in surveying the emission of HI in 
the local Uni verse, such as the HI Parkes All Sky Sur- 
vey, HIPASS (JBarnes et al.l I2OOII : iMeyer et all |2004) an d 
the HI Jodrell All Sky Survey, HIJASS (|Lang et all 120031 '). 
The current ly ongoing Arec i bo Le gacy Fast ALFA survey, 
ALFALFA JGiovanelli et al.llJOOsI '). is producing the deep- 
est all sky HI map s of the local Universe ever created (e.g. 
iMartin et ahlbOlOl '). 

The gas that resides inside the galaxy disk forms a 
multiphase medium in which HI gas forms H2 in the dens- 
est regions. The Atacama Large Millimeter Array (ALMA) 
promises to reveal an unprecedented wealth of information 
H2 an important tracer of star formation. Radiation from 
young stars which form in these regions give rise to ionised 
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bubbles of HIT, emit ting intense Lyman- g radiation, notably 
surveyed in SINGS (JMeurer et albooeh . 

In addition to these impressive surveys are equally im- 
pressive plans for next generation radio telescopes. For ex- 
ample ApertiB MeerKAT^ and ASKAlff] are the pathfinder 
missions for the ambitious Square Kilometer Array (SKA) 
which may detect 10^ gala xies out to z ~ 1.5 yielding strong 
cosmological constraints (jAbdalla fc Rawlingg |2005| ). The 
ASKAP al l sky survey WAL LABY, will detect 5 x 10^ galax- 
ies in HI (jDuffv et al.ll201ll ') which will be the largest cat- 
alogue until the completion of the Five hundred metre Ar- 
ray Spherical Telescope (FAST) which complements SKA in 
the Northern hemisphere, po tentially detecting 10^ galaxies 
with mean redshift z ^ 0.15 JDuffv et aLlbOOar ). 

Together these planned telescopes will revolutionise our 
ability to measure both the local HI mass function and its 
evolution to high redshifts, as well as the distribution of 
HI within galaxies across a large fraction of the Universe. 
The expected HI data is, however, far in advance of current 
theoretical models due to the intrinsic complications of HI 
formation and destruction in galaxies. This is in addition 
to the already uncertain gas treatment in the theoretical 
models for sub-grid baryonic physics. 

Although the detailed distribution of HI is not well un- 
derstood, the basic picture of the lifecycle of HI can be rea- 
sonably well described through simple physical arguments. 
The gravitational potential of the dark matter forms the 
framework in which gas can accumulate, cooling through a 
variety of thermal and non-thermal emission mechanisms. 
However, due to the cosmic UV background only haloes 
with a virial temperature greater than ~ 10* K can retain 
this gas, thereby setting an effective l ower mass floor fo r 
the haloes that can host galaxies (e.g. iQuinn et al.lll996r ). 
Eventually, as the gas density increases, it can become self- 
shielded against this background UV flux and will form large 
collections of cold, neutral hydrogen. In the centres of these 
clouds molecular hydrogen may form which will lead to a 
rapid cooling of the gas. The dense core of clouds may then 
experience therm o- gravitationa l instabilities which can lead 
to star formation l|Schavell2004l ). The high-mass stars will in 
turn ionise the very regions that helped create them. This 
is accomplished directly via photo-ionising radiation from 
the stars, as well as in driving bulk motions of the gas from 
the galaxy, in turn reducing t he self-shie lding against the 
background radiation field fe.g. IPawlik fc S cliavc 2 003 ) . 

In an ideal simulation one would trace the impact of 
individual stars, or stellar regions, in ionising the gas and 
thereby accurately track the non-local effects of star forma- 
tion on the Interstellar Medium (ISM). This is an extremely 
computationally expensive process for the large number of 
galaxies formed in hydrodynamical simulations of cosmolog- 
ical volumes, such as is considered in this work. However, 
by limiting the radi ative transfe r sche me to postprocess- 
ing existing datasets, lAltav et al.l (|20ld ) were able to accu- 
rately model the neutral hydrogen column density distribu- 
tion fu nction at z — 3 {a, separate study bv lMcQuinn et al.l 
(|201ll ) was similarly successful using a postprocessing ra- 



diative transfer scheme at z = 3). lAltav et al.l (|20ld ) found 
that it was critical to model self-shielding above column den- 
sities of 10^* cm~^ and molecular hydrogen formation above 
lO'^^'^ cm~^. In particular we adopt their empirical treat- 
ment of molecular hydrogen formation at high gas densities 
as discussed in Section [3.21 

Several other works have stud i ed the distribution of HI 
at z = 3 such as Pontzen et al.l (|2008l ). iRazoumov et al.l 
(|2008l ) and iTescari et al.l ( 2003 ) but most r elevant for our 
investigation of the low redshift Universe is IPopping et al.l 
(|200g ). henceforth P09. We argue that our higher resolu- 
tion simulations allow us to make more realistic and better 
resolved predictions for the HI, H2 and cold gas mass func- 
tiouQ extending this work to the evolution of the mass func- 
tions and investigate the sensitivity to the physics scheme 
incorporated; a significant improvement on all existing hy- 
drodynamical studies to date. 

Other works investigating HI in galaxies across cosmo- 
logical s cales have taken a semi-a nalytic approach, for ex- 
ample in lObreschkow et al.l ((20091) they calculated that the 
neutral hydrogen in galaxies was distributed across atomic 
and molecular states with a ratio, R — A/hj/A/hi that in- 
creased with the gas mass of the halo. Adopting this ap- 
proach they found that semi-analytic models could repro- 
duce the z = HI mass function above Mhi ~ 10^ M0 , 
below this mass they found a factor two over abundance of 
galaxies relative to the observations. This was opposite to 
another semi-anal ytic study b a sed o n the same underlying 
DM simulation; in I Power et al.l (|20ld ) it was noted that such 
low mass HI s ystems were t oo rar e. Interestingly it has been 
suggested by iLagos et al.l (|201ll ) that if one were to split 
the cold gas into HI and H2 during the simulation, rather 
than in postprocessing, and form stars based on the latter 
alone, then many more low mass HI systems survive in the 
simulati ons. 

In lObreschkow fc Rawlingj ((20093) the ratio R 



from [Obreschkow et al.l ( 20091 ) was applied at higher red- 
shifts. With this scheme they found that there was lit- 
tle evolution in the cosmic density of HI in the Universe, 
JIhi, to z = 1 beyond which the value dramatically de- 
creased, with only half of the local HI predicted to exist 
in galaxies by 2 = 3. This is in contrast to conclusions made 
on the basis of observations of dampe d Lya systems (e.g. 
[Prochaska et al.|[2005l : [Rao et al.|[2006l ) which infer a factor 
two i ncrease by z ~ 3. Further observations from DLA sys- 
tems ((Prochaska fc Wolfe|[2009l ) indicate that the comoving 
HI density at 2 = 2 is, within the errors, co nsistent with 
the r esult at z = as measured by HIPASS ((Zwaan et al.[ 
(20051) (with a best-fit value only 5% higher). This local 
value is similar to the partly completed ALFALFA sur- 
vey of the local Universe, although they find a value 16% 
higher ((Martin et al.[[2010l '). Tantalising evidence that the 
HI cosmic density at z = 0.8 is also similar to that at z = 
is given bv (Chang et al.( ((2010(), who cross-correlated galax- 
ies from optical catalogues with brightness temperature ra- 
dio maps to measure a st atistical 'stacked' es t imate of the 
HI emission. Furthermore. (Prochaska fc Wolfg ((2009|) found 
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* Notably, we find that their mass function studies consider mass 
ranges we find to be unresolved, we demonstrate this in Fig. IA3I 
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that this value increases weakly beyond z = 1, only doubling 
over the redshift range z ~ 2 — 4. Therefore, although the 
integrated HI mass function is observed to be roughly con- 
stant in time, the distribution of HI in mass may vary. This 
is the main subject of our work and the goal of upcoming 
deep radio surveys. 

This paper is organised as follows. In Section [5] we 
describe the physics models and range of resolutions that 
we have simulated as well as the numerical tests we have 
performed to remove artefacts. In Section |3] we detail our 
methodology of converting the hydrogen abundances of the 
simulation into HII, HI, H2. We analyse the HI properties 
of haloes in Section 3] First we consider the distribution of 
the stellar mass and gaseous components as a function of 
halo mass in Section [4.21 We then probe the HI and H2 as 
a function of stellar mass in Section 14.31 As a last probe of 
the distribution of the baryons within haloes we examine 
their half-mass radii for various halo masses and redshifts in 
Section [4.41 We then consider the HI and H2 mass functions 
and their evolution with redshift in Section [S] In Section [6] 
we vary the sub-grid physics schemes to probe the sensitiv- 
ity of the mass function predictions to the particular physics 
prescription. Finally, we conclude in Section [T] 



2 SIMULATION DETAILS 

We have made use of a number of high resolution, 
large-volume hydrodynamical runs created as part of 
the Overwhelmin gly Large Simulations project (OWLS; 
ISchave et al.ll2010l ). This series of simulations all share the 
same initial conditions, i.e. are identical simulation volumes, 
but are resimulated using ~ 50 different physics models 
which allows robust testing of the relative importance of 
each scheme. OWLS has seen a number of successes in repro- 
ducing observed galaxy pro perties. For example , the cosmic 
star formation hist ory (e.g. Madau et al.lll996|') was shown 
to be repro duced in Schave et al. J2010|) . In lBooth fc Schavd 
l|2009l ) and iBooth fc Schavd (|201of rthe AGN implementa- 



tion was found to recover a number of empirical proper- 
ties of super-massive blac k holes, such as the black hole 
mass - halo mass relation (JBandara et al.l 120091): the blac k 
hole mass - stellar mass relation (e.g. iHaring fc Rixll2004h : 
and the black hole mass - bulg e velocity dispersion correla- 
tion fe.g. lTremaine et aLll2002l 'l. Furthermore, this AGN im- 
plementation is capable of reproducing a number of galaxy 
group properties such as their X-ray temperature and den- 
sity pro files as well as galaxy s tellar masses and age distri- 
butions (|McCarthv et al. 2010l). O f partic ular note for this 
work is the study by lAltav et al.l (|201(J ) who found that 
OWLS can reproduce the observed z = 3 HI column density 
distribution function, when post-processed with a radiative 
transfer model. More generally, a particular strength of the 
OWLS suite is that the variety of models allows us to test 
the sensitivity of a particular observable phenomenon to the 
details of the physical processes included (e.g. metal enrich- 
ment, feedback mechanisms etc). 

The OWLS simulations were run with an up- 
dated version of the pu blicly- available GADGET-2 A''- 
body/hydrodynamics code ( Springe] 20051') wi th new mod- 
ules f or radiative cooling (IWiersma et al.ll2009r ). star forma- 
tion (jSchave fc Dalla Vecchial I2OO8I ) , stellar evolution and 



mass loss (IWiersma et al.l 20091). ga lactic winds driven by 
SNe llDaUa Vecchia fc Schavd 1200^ 1 and, for one model, 
AGN (JBooth fc Schavdl2009i r^Each simulation employed A^^ 
dark matter (DM) particles, and N^ gas particles, where 
N — 128, 256 and 512, within cubic volumes of comoving 
length, L = 25, 50 and 100/i~^Mpc. For all cases, glass-like 
cosmological initial conditions were generated at z = 127 us- 
ing the Zeldovich approximati on and a transfe r function gen- 
erated using CMBFAST (v. 4.1. ISeliak fc Zaldarriaga- 199QV 

Simulations with the smallest box were run to z = 2 
and the two larger volumes were simulated down to z = 0. 
For the present study the complete range of box sizes exist 
fo r only one of the p hysics models, that termed 'ZC-WFB^ 
in iDuffv et al.l (|2010f ) . Table [l] lists the available box sizes 
and other pertinent numerical parameters, including the 
Plummer- equivalent softening lengths and particle masses. 
(Note that the dark matter particles have a mass that is 
a factor (1 — /b"'^)//)^"'^ greater than the mass of the gas 
particles, where /""'^ = flh/^m is the Universal baryon frac- 
tion.) 

We use th e Wilkinson Microw ave Anisotropy Probe 
3 year results (jSpergel et al.l 120071 ) , henceforth known as 
WMAP3, to set the cosmological parameters [Qm, Ob, JIa, 
h, 0-8, ris] to [0.238, 0.0418, 0.762, 0.73, 0.74, 0.95] and 
juniv _ QiYQ. Although we use a cosmology that has a 
(78 a bout 10% lower tha n the current year 7 WMAP re- 
sult (JKomatsu et al.ll201lh we believe that the effects on our 
predictions for the various phases of hydrogen will be small. 
This is because only the high-mass end {M ^ 10^* /I'^M©) 
of the mass function is exponentially sensitive to a change 
in as whereas the majority of the HI signal in our simu- 
lations is due to lower mass DM haloes near the 'knee' at 
M ~ 10^ h~^MQ which are much less sensitive to as- When 
comparing DM only simulations with ag values for WMAP3 
and 7, the overall number counts agree within 5 per cent 
for haloes of mass M = 10^^ - 10^^ /I'^M©, while clusters 
{M = 10" - 10^^ /i"^Mq) are nearly 70 per cent more com- 
mon in the higher as universe. We note that the study of the 
dependence of the hydrogen phases on the sub-grid physics 
will be insensitive to the cosmology used, as all simulations 
have the same initial conditions. 



2.1 Simulation models 

We briefly recap the phys ics prescriptions use d in this work 
but direct the reader to ISchave et al.l l|2010l ) and the ref- 
erences therein for a more extensive explanation. We are 
primarily interested in the sensitivity of the HI abundance 
to differing feedback schemes and identify such simulations 
with the labels 'WFB\ 'SFB' or 'NFB' to denote stellar 
feedback that is weak, strong or absent, respectively. Fur- 
thermore we denote a simulation that includes AGN feed- 
back by 'AGN'. Lastly, we want to test the sensitivity of 
the results to the additional cooling of the gas due to the 
presence of metals and denote these simulations with 'ZC\ 
while those without are labeled 'PrimC. Refer to Table [2] 
for an overview. 

For the simulations labeled 'PrimC we model the cool- 
ing of gas via line emission from primordial elements, H 
and He. Those labeled 'ZC have additional contributions 
from nine metals: C, N, O, Ne, Mg, Si, S, Ca and Fe. In 
both cases, additional cooling via free-free Bremsstrahlung 
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Table 1. Details of the various simulation parameters used in this study. From left to right the columns show: simulation identifier; 
the physics scheme(s) available; comoving box-size; number of dark matter particles (there are equally many baryonic particles); initial 
baryonic particle mass; dark matter particle mass; comoving (Plummcr-c quivalent) gravitatio nal softening; maximum physical softening; 
final redshift. Note that we have adapted the OWLS simulation list from lSchave et al.1 1120101 ) such that REF in that work is referred to 
as ZC-WFB here, while All corresponds to the full range of simulations listed in Table |2] 
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emission is also present. The cooling r ates for the different 
elements are taken from 3-D CLOUDY (|Ferland et al.lll998l ) 
tables which assume collisional equilibrium before reionisa- 
tion [z > 9) and photoionisaton equilibrium afterwards, due 
to the p resence of an evolvinR m eta-galactic UV/X-ray back- 
ground (jHaardt fc Madaull200lh . For more details of the im- 
plementatio n of gas cooling (and h eating) rates in the sim- 
ulations, see lWiersma et al.l (|2009r ). 

Gas particles are converted stochastically into stellar 
particles (such that particle numbers are conserved) at a 
pressure-dependent rate in acco rdance with the Kennicutt- 
Schmidt law (f or a review see iKennicuttI 119981 ) using the 
prescriptions of lSchave fc Dalla Vecchial (I20p8|). We assume 
a Chabrier initial mass function (IMF; see IChabrieij |2003| ) 
within each star particle and track the met al mass released 
by the stars using the method described in IWiersma et al.l 
(l2009l) . 

The stars formed in the simulations inject energy in 
neighbouring gas particles via kinetic feedback. For our as- 
sumed IMF, the energy per solar mass of material available 
from supernovae (SNe) is 1.8 x lO^'^ergMQ^. This is dis- 
tributed amongst the local gas according to the dimension- 
less mass-loading parameter 77, in which a total gas mass of 
r/MsN is given an initial wind velocity of v„. This results in 
a fraction 



u 



:0.4(? 



600 km 6 



esN 



1.8 x 10«ergM" 



(1) 



of the energy of the SNe being imparted in the form of 
kinetic feedback (|Dalla Vecchia fc Schavell2008l ): the values 
r) = 2 and Ww = 600kms~^ were adopted for the ^WFB'' 
model. For the ^SFB' simulations Ww {ij) increases (de- 
creases) with the gas density in the region of the star while 
keeping f^ fixed. Note that the kinetically-heated gas parti- 
cles are never decoupled from the hydrodynamics and there- 
fore will tend to collide and 'drag' a larger number of gas 
particles with them, resulting in an effective mass loading 
which is higher than implied by 77. 

In the simulations labeled 'AGN' the black hole at the 
centre of a halo grows at the expense of the surrounding 
gaseous medium by gas accretion as well as through merg- 
ers between black holes. The maximum accretion rate is 
capped at the Eddington limit and 1.5% of the rest mass 
energy of the accreted matter is distribu ted thermally to 
the nearby gas particles as explained in iBooth fc Schavj 



(|2009l ). iBooth fc Schavd (|2009l ) demonstrated that for this 
value of the feedback efficiency the observed scaling relations 
for supermassive black holes are reproduced. 



2.2 Identifying Galaxies and Haloes 

In this work we will consider the HI content of galaxies and 
the associated HI gas lying between them, which is par- 
ticularly important in groups and clusters. We initially de- 
termine the halo in a simulation u sing a Friends-of-Friends 
(FoF) algorithm (jDavis et al.ll 19851 ) with a linking length b. 
This will link all DM particles together that are within a 
distance bl, where I is the mean inter-particle distance in 
the simulation. The gas is then attached to the nearest DM 
particle. The halo can thus have arbitrary shape, defined 
by an iso-density curve of overdensity relativ e to the back- 
ground p/p « 3/(27r63) (|Cole fc Lacevlll996l ). We choose a 
value b = 0.2 which, for an isothermal densi ty profile, will 
have a halo mean over-density of {p)/p ~ 180 (jLacev fc Cold 
[l99J). This is close to the spherical top-hat collapse model 
prediction for a virialised object « ISvr^ ~ 178. 

In the simulations with L = 100/i~^Mpc we typically 
resolve 230 (15) haloes of virial mass lO" - lO" /j-^Mq (> 
10" h-^Me) a.t z = and for L = 25 h'^Mpc there are 840 
(13) haloes of 10^^ - lO" /j-^Mq (> lO" /i-^Mq) a.t z = 2. 
As all simulations of the same volume share the same initial 
conditions, the haloes formed are identical, except for the 
specific choice of sub-grid baryonic physics. This means that 
the number of haloes between simulations typically differs 
by less than 5 per cent, essentially due to slight variations in 
mass of specific haloes which can affect whether they meet 
mass cuts or not. 

To identify gravit ationally bound s tructures within 
haloes we use SUBFIND (|Dolag et al.ll2009l ). This allows us, 
for example, to identify gas that has been stripped from an 
infalling galaxy but which is still bound to the background 
cluster, as well as satellite galaxies. Unless stated otherwise, 
the masses and sizes which we will quote are for objects 
found with subfind (and we will refer to such objects as 
galaxies), but when we discuss the halo in full (including 
satellites) we will use that which comes from the FoF algo- 
rithm and will add the 'Halo' superscript to the variable. 
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Table 2. A list of the simulated sub-grid baryon p hysics scheme s utilis ed in this study. From left to right, the columns list the simulation 
name used in this study; the name as defined bv lSchave et alj pOlOh : and a brief the description of the physics modelled. The list of 
available simulations run with each model (box-size and resolution) is given in Table [l] 



Simulation 



OWLS name 



Brief description 



PrimC-NFB NOSN^NOZCOOL No energy feedback from SNe and cooling assumes primordial abundances 

PrimC-WFB NOZCOOL Cooling as PrimC-NFB but with fixed mass loading of winds from SNc feedback 

ZC-WFB REF SNe feedback as PrimC-WFB but now cooling rates include metal-line emission 

ZC.SFB WDENS As ZG.WFB but wind mass loading and velocity depend on gas density 



ZC.WFB.AGN AGN 



As ZG.WFB but with AGN feedback 
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Figure 1. A simple flow diagram indicating the methodology to 
attribute the states of hydrogen to the gas particles of the simu- 
lation. In our model, the HII and H2 components arise from the 
diffuse (intergalactic) and star-forming (interstellar) components 
respectively, whereas the HI signal is found in all regimes. Note 
that we only ever consider particles that belong to gravitationally 
bound structure. 



3 MODELLING THE NEUTRAL HYDROGEN 

In this section we discuss our method for calculating the 
mass fraction of hydrogen in neutral molecular (H2), atomic 
(HI) and ionised (HII) phases, for each gas particle. The to- 
tal hydrogen mass fraction is first estimated for each particle 
by taking a smoothed average over its neighbours, using the 
SPH kernel (a value for the hydrogen abundance is carried 
by individual particles). The methodology for splitting the 
gas into the 3 phases of hydrogen depends on whether it 
is low-density intergalactic gas or high-density interstellar 
gas. We consider each one in turn below, with the overall 
procedure summarised as a flow diagram in Fig. [T] 



3.1 Intergalactic Gas 

Most of the gas resides outside, or in the outskirts, 
of galaxies and will have a significant HII compo- 
nent, due to both the photo-ionising UV/X-ray back- 
ground and coUisional ionisation resulting from gravita- 
tional shock heating. OWLS tracks an evolving UV/X-ray 
background (jHaardt fc Madaul |200J ) after reionisation at 
z = 9, that primarily affects the cold, low-density inter- 



galactic medium. Gas at higher density can cool more effi- 
ciently, leading to larger neu tral (HI) fractions (for full de- 
tails see lWiersma et al.ll2009^ . A significant fraction of gas is 
shock-heated to much higher temperatures at low-redshift, 
creating the warm-hot intergalactic medium and the intra- 
cluster medium. We use pre-calculated i onisa tion tables for 
hydrogen from CLOUDY JFerland et al.l 11998 '). to compute 
the optically thin HI and HII fractions of gas particles as a 
function of their current temperature and (SPH smoothed) 
density. 



3.1.1 Optically Thick Gas 

For gas that resides at moderate densities (in the range 
TiH ~ 10~* — 10^^ cm^^) assumptions of optically thin gas 
exposed to a uniform UV background start to break down. 
Firstly, the local galactic UV intensity can be significantly 
stronger than the global average, from the presence of young 
massive stars or a quasar (the so-called proximity effect). In 
this case, the CLOUDY tables will overestimate the amount 
of HI present as the assumed photo-ionisation rate will be 
lower than it should be. Secondly, near galaxies the gas den- 
sity can reach levels which are sufficient to self-shield the 
gas from the global UV background, which in turn would 
lead us to underestimate the fraction of HI. As we have pre- 
viously discussed, ideally we would perform a full radiative 
transfer calculation to distribute ionising flux from nearby 
stellar sources through the gas. Such a technique would also 
lead to a more realistic calculation of the self-shielding in 
dense regions. We turn to simpler prescriptions that can 
be more easily applied to large volume simulations. In this 
study, we assume that only the self-shielding is important 
and attempt to correct for it using a simple method guided 
by the observational data. 



3.1.2 Self-shielding Correction 

We follow a methodology similar to that of P09 and choose 
to approximate the onset of self-shielding by adopting a 
pressure-based threshold. It is reasonable to expect that self- 
shi elding occurs above a given surf ace density. As, shown by 
e.g. ISchave fc Dalla Vecchial (|2008n , the pressure in the disk 
ought to scale quadratically with the surface density as the 
thickness of the disk is set by the local Jeans length. We 
can even use similar arguments to deflne an order of mag- 
nitude estimate for the pressure threshold value PshieW by 
considering the local volume density at the threshold, riahieid, 
related to the pressure by PsMcid/k = (Tishieid/fc)!^ K with 
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Table 3. A list of the pressure thresholds used to approximate 
the onset of self-shielding, Pshiold/fci for gas particles that have 
temperatures below 10^'^ K in methods SSalfalfa (2nd column, 
normalised to low-redshift ALFALFA observations) and SSdla 
(3rd column, normalised to high-redshift Damped Lyman-a ob- 
servations). Values in brackets denote pressure thresholds that 
produce HI densities 10 per cent higher and lower than observed, 
respectively. Zero values correspond to cases where even mak- 
ing all cold gas in the galaxy neutral gives insufficient HI to 
match the observations (typically ~ 20 per cent below the de- 
sired value). The neutral interstellar gas is still separated into 
atomic and molecular neutral components. 



Simulation 



SS^ 



SSr 



PrimCNFB 

PrimC. WFB 

ZC.WFB 

ZC.SFB 

ZC.WFB.AGN 



[0, 0] 

[0, 0] 

150 [133,175] 

74 [60, 93] 

27 [21, 35] 



159 [132, 193] 

47 [30, 68] 

98 [80, 123] 

17 [3, 28] 

76 [62, 93] 



a calculation of nshidd from Equation 11 in ISchayg (|2001bf ) 



"-H, shield 



l + z\-9/Vnb/l'^"'^' / f- ^-'/' 



0.02 



0.16 



iVni 



2.7 X 10" ern- 



es) 



giving riH, shield — 10"'^ — lO"'^ cm"'^ for 2 = 3 and assum- 
ing that overdensities relative to the background density, S, 
of order 10^ represent self-shielded galaxies, which have a 
fraction /g of their disk mass in HI. This is equivalent to 
a pressure Pshieid/fc ~ 10 — lO^Kcm^'^. This simple esti- 
mate h as been con fi rmed using a radiative transfer simula- 
tion in lAhav et"aLl (|201G| ). 

The pressure threshold used in P09 was tuned to re- 
produce the mean HI den sity from the 100 Brightest 
HIPASS galaxy catalogue (jZwaan et al.l 120031 ). finding a 
value Pshieid/fc = ISSKcHi^^. They only applied this to gas 
cool (or dense) enough to have a recombination time less 
than or equal to the sound crossing time (see Fig. 2 in P09). 
While it is clear that hot, coUisionally ionised gas should be 
excluded, the physical motivation for the particular criterion 
chosen by P09 is unclear to us. We simply use an upper tem- 
perature limit for gas that is considered to be self-shielded, 
essentially forming a 'cold' gas population with T ^ 10^'^ K. 
Our a pproach is similar to that adopted by I Power et al.l 
(|2010f ) for their investigation of cold (T ^ 10^ K) gas mass 
functions from semi-analytic simulations. 

In what follows we have treated the self-shielding regime 
in three different ways. In what we shall call SSNone we 
shall presume that there is no self-shielding. This is clearly 
unphysical but it can be helpful to compare this with 
the other two methods. For method SSalfalfa we match 
the measured cosmic HI density in the HI mass range. 
Mm = 10^° - W^^h-^M p,, from the 2 = ALFALFA sur- 
vey (JMartin et al.ll2010^ . This normalisation is computed 
for the 100 /i~^Mpc simulation as this volume has been run 
with all sub-grid physics models. Cold gas above the pressure 
limit is then considered to be completely neutral, while gas 
below the threshold has HI and HII fractions taken from 
the optically-thin CLOUDY tables discussed previously. We 



find that a pressure of Pshieid/fc ~ 150 Kcm ^ matches the 
ALFALFA result for our defauh model {ZC.WFB, which is 
also the simulation closest to that studied by P09 in terms 
of the sub-grid physics modelled). The pressure thresholds 
used for all simulation models are listed in Table [3] typically 
in runs with stronger feedback we require a lower pressure 
threshold to match the HI density while in the runs with 
primordial cooling and weak or no feedback, there is insuffi- 
cient (~ 20 per cent too low) HI gas in simulated galaxies to 
match the observations at low-redshift (the reverse is true 
at high-redshift). 

There is no reason to believe that the pressure thresh- 
old will be independent of redshift. In fact, we expect 
it to de pend on the i ntensity of the evolving UV back- 
ground (|Schavd l2001bl ). For method SSdla we therefore 
normalise Pshieid to agree with observations at higher red- 
shift, namely the cosmic HI density inferred from Damped 
Lyman-a absorption sy stems. In particular, we wi ll use mea- 
surements presented in lProchaska fc Wolfd ([20091) for 2 « 2. 
These results are close to the 2 = HI cosmic density in- 
ferred by ALFALFA (e.g. Hhi.dla = 1.05r2Hi,ALFALFA = 
3.75 X 10~*), but lead to somewhat higher values for Pshieid. 
We tune the pressure value until the integrated HI mass 
function (where available, as shown in Fig. Ill|l matches 

f2HI,DLA. 



3.2 Star- forming Interstellar gas 

As gas cools, its density will eventually increase until its 
structure can no longer be modelled reliably in our sim- 
ulations. Such gas, associated with the star-forming inter- 
stellar medium in galaxies, is expected to be multiphase 
for hydrogen particle num ber densities above n^ ~ lO"'^ 
- IQ-^an-^ (|Schave|[200i ). In OWLS, such gas particles 
are identified when they exceed a critical particle density, 
njj = 0.1 cm""^ and have a temperature smaller than 10^ K, 
after which the gas density and pressure is evolved according 
to an effective polytropic Equation of state 



P — Pcrit I — — 

n„ 



(3) 



where Pcrit/fc = 2.3 x lO^Kcm"^. All runs considered as- 
sume 7eff = 4/3 which has the a dvantage that the Jeans 
mass is independent of density (jSchave fc Dalla Vecchial 
[2OO3), avoiding the situation where particle masses are 
greater than the Jeans mass in dense regions, a numerical 
artefact k nown to result in arti ficial fragmentation of the 
gas cloud (JBate fc Burkertlll997^ . 

The Equation of state gas, which we will refer to as in- 
terstellar gas, lies in dense regions of the simulation where 
self-shielding effects are important. We therefore assume 
that the hydrogen exists as either HI or H2 (we ignore in- 
terstellar HII regions as they will contribute little to the 
hydrogen mass in galaxies). To separate the two phases, we 
utilise the empirical ratio relating H2 and HI surface densi- 
ties, Psurf ~ SHa/Sni, to the local ISM pressu re measured 
in th e THINGS survey (Table 6 and Fig. 17 of lLerov et al] 
l2008h . given by 



Psurf 



P/k 



lO-i-^s Kcm- 



(4) 
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which was observed for a pressure range of 10'' ^ 
P/klKcm"'^] ^ 10^. It is then a trivial step to compute 
the mass of interstellar gas in atomic HI, 



Mm = Mh/(1 + flsurf) . 



(5) 



Note that we have assumed that the HI and H2 in the same 
disk are measured locally. 

In using Equation[S]we must consider whether the pres- 
sure used in the empirical law of Equation |4] corresponds to 
the pressure averaged over the physical length scales that the 
simulations resolve. The THINGS survey reduces its sample 
of 23 dwarf and gas spiral galaxies to a standard resolution of 
~ 20 arc seconds. However, the dwarf galaxies are typically 
closer than the spirals, as they are less luminous. The re- 
solved physical lengthscales for the dwarf and spiral galaxies 
are around 400 and 800 pc respectively. This is close to the 
Plummer-equivalent softening length-scale (in proper units) 
of the TV = 512'' simulations, as given in Table [T] Only the 
100 h^^ Mpc boxsize may be cause for concern when we con- 
vert the pressure to HI mass, but as shown in Appendix lAll 
this simulation appears sufficiently well resolved. 

We also note that Equation [4] is for the local Universe 
only and the assumption that it holds at higher redshifts is 
potentially an i ssue. It is belie ved that this ratio will depend 
on m etallicity (|Schave|[2001al . [2004: Krum holz ct al. 2008, 
2003) which decreases wi t h redshift (e.g. iTremonti et al.l 
20041 : ISavaglio et all l2005l : lErb et all I2OO6I ). We therefore 
caution the reader that the molecular fraction is likely to be 
a more complicated function of galactic properties than as- 
sumed here, given the limited physical information we have 
for galaxies in our simulations. 



3.3 Resultant phase diagram for a specific case 

In Fig. [2] we show the overall neutral fraction in atomic (HI; 
top panel) and molecular (H2; bottom panel) states as well 
as the ionised HII component (middle panel) in the temper- 
ature - density phase diagram at z = and 2 (left and right 
columns respectively) for all gas particles in a simulation. We 
show results from the ZC^WFB 100 h~^Mpc simulation, us- 
ing method SSalfalfa for the self-shielding correction. (The 
equivalent phase diagrams for the other simulation schemes 
are broadly similar.) 

The main result from the top panels is that the majority 
of HI resides near the self-shielding limit, demonstrating the 
importance of such a correction. Furthermore, even though 
at z = the majority of hydrogen resides in the hot plume 
of intergalactic material (middle panels), this gas is typically 
too highly ionised to be a large contributor to the overall HI 
mass. By design, the molecular component exists only in the 
densest material, which obeys the imposed Equation of state 
power law. Differences in the phases between low and high 
redshift are most apparent in the HI distribution below the 
self-shielding limit. At 2: = 2 (right column) the HI content 
of the gas drops more sharply as the more intense global UV 
background ionises more of the gas. 

The relative importance of each phase in gravitationally 
bound structure (and overall in the Universe) is summarised 
in Table 2] as a function of redshift for the 50/i~^Mpc sim- 
ulation with the ZC^WFB physics scheme and SSalfalfa 
self-shielding model. 



Table 4. A list of relative fractions of the hydrogen phases (HII, 
HI and H2) for the simulation model ZC.WFB in the 50/i~lMpc 
simulation volume using self-shielding threshold SSalfalfa- As 
we are unable to resolve the low mass end of the galaxy distribu- 
tion the HI and H2 values are lower estimates. We consider the 
two cases for the hydrogen phases in both gravitationally bound 
structures and overall in the Universe (latter values are in paren- 
theses). 



Redshift 


HII 


HI 


H2 



1 

2 


0.964 (0.993) 
0.862 (0.983) 
0.736 (0.974) 


0.026 (0.005) 
0.085 (0.010) 
0.167 (0.016) 


0.010 (0.002) 
0.053 (0.007) 
0.097 (0.010) 



4 NEUTRAL GAS IN GALAXIES AND 
HALOES 

OWLS provides a wide range in dynamical mass over which 
the distribution of baryons can be measured. Additionally, 
different feedback and cooling schem es can be tested (e.g. 
ISchave et all I2OIOI : iDuffv et al.l I2OIGI '). In considering how 
the HII, HI and H2, as well as the stars formed from the 
dense hydrogen gas, are distributed across the simulation, 
we must separate the intertwined effects of baryon physics 
with mass scales and evolution. The results in this section 
concentrate on one simulation, ZC-WFB, but we will also 
make comparisons with ZC.WFB^AGN (arguably our most 
realistic model) to understand the possible effects of AGN 
feedback on the phases of hydrogen. 

In Fig. |3] we provide a visual overview of the effects of 
halo mass and redshift on the various distributions of the 
hydrogen phases. We compare haloes of virial mass Afvir ~ 
10 , lO'^ and 10 h~ Mq (approximately corresponding 
to a galaxy, group and cluster, respectively), in columns left 
to right respectively. From top to bottom, the rows show 
density maps for the HII, HI, H2 and stellar components. 
The galaxy halo is at z = 2 while the group and cluster 
haloes are at 2 = 0. 

Immediately apparent is the diffuse and filamentary na- 
ture of the HII in comparison with the dumpiness of the H2 
/ HI components, due to the latter two being associated with 
the cold and dense gas. 



4.1 Resolution Effects 

Before we study the neutral gas content of galaxies and 
haloes in detail, we first investigate the effects of resolution 
on our HI and H2 mass determinations. This is presented in 
Appendix|A]where we test the effect of finite mass resolution 
and simulation box size on the predicted mass of HI (Ap- 
pendix |Xl|, the HI mass function (Appendix IA2|) and the 
molecular hydrogen mass function (Appendix I A3[) . In sum- 
mary, we find that the properties of the baryonic species 
in haloes are converged above a certain limiting mass, with 
each species, X, as 



M„ 



= Mlin 



L 



512 



50h-iMpc iV 



(6) 



where A'' is the cube-root of the number of Dark Matter par- 
ticles and L is the comoving length of the simulation cube. 
The co-efficient Miim depends on the species X; in ascend- 
ing mass order these are Mum = [1,5, 5, 10, 500] x 10^ h~^MQ 
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Figure 2. From top to bottom, we present the mass fraction of HI, HII and H2 for all gas in a simulation for a given gas temperature 
and hydrogen particle number density, in the ZC-WFB simulation. Left panels show results at z = and right panels for z = 2. Method 
SSalfalfa is assumed for the self-shielding correction (tuned to the ALFALFA HI mass function at low-rcdshift). Additionally, we have 
projected each fraction as a function of density in the histogram below each image (black line), with the red line a simple cumulative 
distribution function of this histogram. The pressure threshold for the onset of self-shielding is especially clear in the top panels as a 
discontinuity in the HI fraction. The H2 is constrained, by design, to reside in the high-pressure multiphase ISM, modelled as a polytropic 
Equation of state. 



for X= [HI, Cold Gas, Stellar, H2, Virial] respectively. This 
stellar mass limit implies that approximately 500 star par- 
ticles are required for convergence, while the Virial (or To- 
tal) mass limit is dominated by the Dark Matte r component 
and c orresponds to nearly 7000 particles; from iDuffv et al.l 
l|2010l ) we are confident that the masses and profiles of such 
systems are well converged. For HI, this value is more con- 
servative than assumed by P09 who did not present con- 
vergence tests (see Fig. IA3I for the comparison of the P09 
results and our own resolution test). 



4.2 Halo mass fractions 

We now proceed to investigate the mass fractions of the 
various baryonic species in haloes. All available simulation 
volumes listed in Table [l] are used for the two sub-grid 
physics models considered, ZC.WFB and ZC.WFB.AGN. 
Taken together, we are able to probe a virial mass range 
of 5 X 10^° — 5 X 10^* /i~^M0, representing the largest dy- 
namic range in mass yet considered in such a study. Such a 
large dynamic range allows us to probe whether there exists 
a characteristic mass scale at which HI is particularly abun- 
dant. This is an important issue for understanding the selec- 
tion effects that a blind, flux-limited HI sample would suffer 
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Mvir = lO^^/j-lMp, z = 2 



Afyir = 10"fe~^MQ, ^ = 



j\/vii. = 10^'*/i-^Mq, 2 = 



HII 



HI 



Ha 



Stars 




Figure 3. The distribution of various baryonic components in a typical z = 2 galaxy (left), z = group (middle) and z = cluster 
halo (right), with Mvir ~ 10^^, 10^^ and 10^^ /i~^Mq respectively. We show the column densities of different baryonic components in 
the rows; from top to bottom these are HII, HI, H2 and stars. Results are taken from the A^ = 512 simulation with L = 25, 50 and 
100/i~^Mpc, left to right columns respectively. Note that different panels show different size regions, as indicated in each panel. The 
coordinates (column densities) are in comoving (proper) units. Red and blue circles denotes -Rvir and O.lRvin respectively. The boxes 
have the same length along the line of sight as across, as indicated by the arrow. All images have been rotated to present the object 
within O.liJvir to be face-on. 
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from, as well as the prospects of using HI surveys for cos- 
mology. We only consider self-shielding method SSalfalfa 
in this section as we are primarily interested in the sensi- 
tivity of the various components to halo mass, redshift and 
galaxy formation modelling (specifically whether we include 
AGN or not). 

In the top row of Fig. [4] we present the mass frac- 
tion of the baryonic species for ZC-WFB, as a function 
of halo mass at redshifts z = 0,1 and 2 (left, middle and 
right panels, respectively). The bottom row shows the corre- 
sponding results for the simulation including AGN feedback 
{ZC.WFB_AGN). 

At z = the ionised gas is the largest baryonic mass 
component of all haloes, followed by stars, although the two 
are comparable for Afvir ~ 10^^ - lO" H^'^Mq. This pattern 
is seen at all redshifts both with and without AGN, although 
the ionised component is less dominant at higher redshifts. 

For the simulation without AGN feedback (top row) the 
HI has a more complex distribution. At z = 0, top left panel, 
there is an initial rise and then a clear decrease in the HI as 
a fraction of the total halo mass, beginning at 0.5 percent at 
M^if° ~ 10^^ h~^MQ and decreasing to nearly an order of 
magnitude less at M^,f° « 5 x 10^"* /I'^M©. The pattern is 
similar at 2 = 2 where the increased dynamic range allows us 
to identify a broader peak around a smaller mass, M^^'° ~ 
5x10^^ h~^MQ. However, at z = 2 the HI mass fractions are 
typically an order of magnitude higher than at z — 0. The H2 
fractions broadly track the HI fractions although they are 
noticeably lower in the lowest mass haloes. In such systems, 
the mass fraction of baryons as a whole decreases as the gas 
is expelled by supernova feedback, leading to less efflcient 
star formation. At higher masses (massive galaxies, groups 
and clusters) the decline in neutral fractions is due to the gas 
being used up and turned into stars, as the supernovae are 
no longer able to blow significant amounts of gas out of the 
halo. Ram pressure stripping of the cold gas in the hot halo 
atmospheres and the larger cooling times also contribute. 

We also see evidence of a form of 'downsizing' (e.g . 
Pozzetti et ai]|2003l : iFontana et"ai]|2004 iDrorv et al1l200l 



2005|) in the HI and H2 fractions as a function of redshift, 
which decrease faster for higher halo masses. We attribute 
this to a runaway conversion of this ga s into stars f or the 
ZC-WFB simulation, as it was noted in iDuffv et al.l (|20ld ) 
that this particular physics scheme overproduces the stellar 
fraction in groups. This can be seen in the bottom panels 
of Fig. |4] where we include AGN which prevents the cooling 
flows, giving a low stellar (and cold gas) mass fraction. AGN 
feedback strongly suppresses the molecular mass fractions. 
The HI fraction, however, is typically slightly greater than 
in the simulation without AGN feedback. This is due to the 
cold gas typically being at lower density in the AGN model 
due to the stronger feedback. 

4.3 Stars and Neutral Hydrogen 

In Fig. [5] we show the HI mass as a function of the stellar 
mass of the galaxies in ZC^WFB at z = 0, 1 and 2, top to 
bottom panels respectively. There is a clear positive correla- 
tion between the two parameters, namely that larger stellar 
niciss systems also have more HI, as shown in the left column. 
We probe this correlation over approximately 2 decades of 
stellar mass at z = 2. 



The correlation seen in the left column of Fig. [5] is less 
than linear however, as is clear from the middle column 
where the ratio of HI to stellar mass is presented as a func- 
tion of the stellar mass. The ratio is a decreasing function of 
stellar mass, which means that, for their size, massive galax- 
ies are relatively deficient in HI when compared with smaller 
galaxies. This is in very good agr e emen t with observations 
at low-redshift from lFabello et al.l (|201ll ). who combined HI 
and stellar data from the ALFALA and SDSS surveys, re- 
spectively (their results are plotted in Fig. [5] as red circles). 
We find that the trend continues at higher masses with the 
final data point indicating perhaps a turnover in the power 
law, a regime which is, however, influenced by the lack of 
AGN in this model. However, as seen in the right column of 
Fig. [31 model AGN typically has more HI at a given stellar 
mass than observed. 

The ratio of HI to stellar mass increases with redshift. 
For example, at z = 2 a galaxy of stellar mass 10^^ h~^yiQ 
will have nearly an order of magnitude more HI than a sim- 
ilar galaxy at the present day. Such an evolutionary trend 
is expected as the accretion of cold gas on to galaxies slows 
down with time and is no lo nger sufficient to replace th e gas 
used in star formation (e.g. Ivan de Voort et al.ll2011al lbr). 

In Fig. [S] we consider the case of the molecular hydrogen 
as a function of the stellar mass of objects, for z = — 2 
for a simulation with SNe feedback, and one with addi- 
tional AGN feedback {ZC.WFB and ZC.WFB.AGN). For 
the f ormer simulation a t z = we compare to local galaxies 
from lLerov et al.l (|2008l ) we find that the simulations match 
the low stellar mass end well over the stellar mass range of 
the observations. 

The situation at 2; = 2 is more complex. At the low 
m ass end the ZCWFB simulation matches the observations 
of iTacconi et al.l ((20101) , however for systems with stellar 
mass ^ 10^^ /i~^Mq the observed galaxies typically have 0.5 
dex greater molecular mass than the simulations. This situa- 
tion becomes even more extreme at z = 1 with the disparity 
at the high stellar mass end now an order of magnitude. We 
can attempt to match these observations by suppressing the 
conversion of cold gas into stars through additional feedback 
with an AGN but the discrepancy increases, with the molec- 
ular mass of simulated galaxies lying below the observations 
at all redshifts. If we compare this with the over-prediction 
of the HI mass for a given stellar mass (bottom right column 
of Fig. [SJ we see that for a given stellar mass, the simulation 
with AGN feedback is probing haloes with a much greater 
halo mass and hence will have more shock-heated gas. In our 
formalism this will move gas from a dense regime in which 
it is dominated by molecular hydrogen into a lower density, 
hotter gas that has more atomic hydrogen, this redistribu- 
tion of gas to larger radii is demonstrated in Fig. [T] 

4.4 Spatial distributions 

As well as the total mass of each species in a given halo, 
wc can probe the spatial distribution of the gas and stars 
within the haloes, as shown in Fig. O To quantify this Fig. [7] 
shows the half- mass radius for each species, Rf/2, plotted as 
a function of halo mass for the different redshifts. Note that 
we only include particles that are gravitationally bound to 
the halo. Again, results are shown for the ZC^WFB runs. 
In haloes more massive than ~ 10^"^ /i~^Mo, the HI 
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Figure 4. The HI, HII, H2 and stellar mass fractions of haloes as a function of the total halo mass for 2: = 0, 1 and 2 , left to right 
respectively. The horizontal dotted line is 0.9 /||^"'^, the value found in non-radiative gas simulations of lCrain et al.l 1(2003) ■ We consider 
runs with metal line emission cooling and supernovae feedback, ZC-WFB, in the top row and then the runs with the AGN scheme added, 
ZC-WFB^AGN, in the bottom row. This latter case is particularly important when considering the stellar fractions in massive haloes. 
Method SSalfalpa (using ALFALFA data) is used for the self-shielding approximation. We only consider haloes with a virial mass 
above the limit in Equation [B] these haloes typically have 10"* dark matter particles. Note the smaller dynamic range in the AGN case 
at z = and 1 where we do not have an L050N12 simulation to probe galaxy scales. 



distribution rapidly becomes more extended than both the 
stellar and molecular H2 distributions, growing to nearly half 
the size of the ionised diffuse component for haloes of mass 
10^"* /i~^Mq. Since the molecular gas is only identified with 
the high-density star-forming gas, it is not surprising that it 
has a similar half-mass radius as the stars (large differences 
start to appear in low-redshift groups and clusters due to 
the presence of a diffuse stellar component). 



5 NEUTRAL GAS MASS FUNCTIONS 

A key observable for each state of hydrogen is its mass func- 
tion. In this section we focus on the HI mass function (Sec- 
tion I5.1|l and the molecular and cold gas mass functions 
(Section l5.2p . for the ZC.WFB simulations as they span the 
largest range in mass. Results are presented at ^ = 0, 1 and 
2 for the three different self-shielding methods. The effects 
of feedback and metal- line cooling will be investigated in the 
following section. 



The power- law dependence of size on halo mass for the 
io nised and stellar components is a generic prediction given 
bv lMo et al.l l|l998l ). in which they also predicted that haloes 
of a given mass will have smaller half-mass radii at higher 
redshift (due to the higher density and hence smaller virial 
radius). We present simple power-law fits to our simulated 
data in Table (5] The half- mass radius of the ionised hydro- 
gen has a cube-root dependency on the halo mass which 
is a generic prediction for a component with a smoothly 
distributed isothermal profile. The molecular hydrogen and 
stellar components are more centrally concentrated than the 
HI and HII at all masses and redshifts. They are sensitive 
to the implementation of the AGN feedback which dramati- 
cally increases the typical half-mass radius by nearly an or- 
der of magnitude. However, at ^ = the effects of the AGN 
feedback in increasing the half-mass radius of the stars and 
H2 is lessened, although the stellar component is, roughly, 
a factor 5 more extended due to the suppression of cooling 
flows which would lead to a highly concentrated, massive 
central g alaxy (a similar effect was noted for galaxies at 
2 = 2 bv lSales et al.ll201Cl using the AGN feedback scheme 
discussed here). 



5.1 HI Mass Functions 

The HI mass functions for ZC_ WFB axe shown in Fig. |8l As 
discussed earlier, we use individual subhaloes to represent 
the galaxies and we determine the minimum HI mass for 
each simulation according to the values given in Equation |6l 
In all cases, the number density of objects decreases as a 
function of increasing HI mass, with a strongly suppressed 
high-mass end and a more gradual slope at the faint end 
(where simulation data is available). The good agreement in 
abundance at the overlap between the different simulation 
volumes, increases our confidence in the numerical resolution 
cuts we have made. 

We consider the no-shielding case first (left column of 
Fig. [8]) and can immediately see that a.t z = Q, top panel, the 
simulation points lie below the ALFALFA results. This justi- 
fies our assumption that the self-shielding correction (which 
increases the amount of HI at moderate densities) must be 
implemented to match observations. In the absence of self- 
shielding there appears to be little evolution in the HI mass 
function between 2 = and z = 2. 

We can tune the self-shielding pressure threshold to 
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Figure 5. In the left column we present HI masses of all subhaloes in logarithmic bins of stellar mass ^ = 0, 1 and 2, top to bottom 
respectively. The middle and right columns show the ratio of HI mass to stellar mass of each halo as a function of stellar mass. Results 
are for the ZCJ\!VFB (left column) and the AGN model (right column, ZCLVFFB.AGN). We have included a range of simulation volumes 
to attain the maximum dynamic range. We find that 5 X 10'^ stellar particles is sufficient to attain a well resolved value. The errors 
are the \u errors estimated by bootstrap analysis and added in quadrature with Poisson error estimates (typically only important for 
the highest mass bins in a simulation). The red circles are from a joint ALFALFA-Sloan Digital Sky Survey comparison at low-redshift, 
compiled bv lFabello et al.l l l201ll 'l. which nicely agree with our ZC.WFB simulation. Intriguingly the more realistic model which includes 
AGN overproduces the observed relation; to attain similar stellar masses the haloes with AGN are more massive. This also means that 
the haloes with have a higher HI mass as shown in Fig. |4] 



Table 5. Power-law fits to the half-mass radius versus halo mass relation for the various baryonic components given by 
logj^Q(i?^, /h^^kpc) = Norm -F Slope X logj^Q(Miiaio/A^pivot) with a pivot mass of 2 X 10^^ /i^^Mq assumed. Quoted errors are 68 
per cent confidence limits, estimated from bootstrap resampling. The top series of values arc for ZC-WFB and the bottom set of three 
redshifts is for ZC.WFB.AGN. 
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HI 

Norm [h.^-'-kpc] 
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Norm [h.^-'^kpc] 
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H2 
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0.52±[i;!J4 
n CI +0.04 

n 4«+0.05 



match the 'knee' of the z = ALFALFA mass function and 
utilise this same pressure cut at higher redshifts, as shown 
in the middle column of Fig.[8l Interestingly, the simulations 
predict relatively little change in the abundance of haloes in 
this mass range with redshift. The high-mass end matches 
the ALFALFA data very well. At the low-mass end, the sim- 
ulation under-predicts the number of objects; this flattening 



is also seen in the cold gas mass function (see Section [5.21) 
and thus appears to be related to the galaxy formation mod- 
elling in the simulation. 



The HI mass functions are well parameterised by a 
Schechter function (particularly at 2 = 2 where we have 
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Figure 6. As in Fig. [5] but now wc consider the molecular mass as a function of stellar mass in the objects for the ZC-WFB (left 
column) and the AGN model (right column, ZC.WFB.AGN) . The errors are the la errors estimated by bootstrap analysis and added 
in quadrature with Poisson err or estimates (typic ally only important for the highest mass bins in a simulation). The re d circles in the 
top ro w are local galaxies from lLerov et al.l 1120081 1 and at higher rcdshifts subsamples (close to the relevant redshift) from lTacconi et al.l 
l|20ld l. 



the highest dynamic range) given by 
dn _ 6 / M 
dM ^ mT 



f-V 



(7) 



with normalisation 6, low-mass slope a and characteristic 
mass scale Mt, respectively. This can be integrated to give 
the mass density 



M 



e A^* dM. 



(8) 



The global HI density (assuming the above masses are HI 
masses) in units of the critical density, pcrit, is then 



^Hi =6»r(a + 2)A/*/pc 



(9) 



where r(x-) is the standard gamma function. Comparing our 
results a.t z — 2 with the low-redshift data, we predict mild 
evolution of the HI mass function with the global HI density 
decreasing from 0.7 to 0.3 at z = 0. This has important 
ramifications for upcoming deep HI surveys such as those to 
be performed with the Square Kilometre Array which will 
detect more HI sources than an extrapolation from the z = 
mass function would suggest. 



Our third method uses a pressure threshold set by 
matching the integrated mass function of the simulations 
at z = 2 t o the DLA inferred cosm ic density of HI at 
z = 2.2 from lProchaska fc Wol5 (|2009l '). The right panels in 
Fig.[S]show that this leads to very similar results as method 
SSalfalfa- The higher global HI density at z — 2 leads to 
a slightly larger value for Mt and a slightly steeper low- 
mass slope, Q. Again, comparing results with the ALFALFA 
data, there is only mild evolution in the global HI density 
with redshift. 

For all models there is however strong evolution seen 
in the low-mass tail, from a steep index of —1.5 at z = 2 
to ^ — 1 at z = 0. This evolution is also seen in the cold 
gas mass function in Fig. [TD] which indicates that it is not 
subject to the self-shielding meth odology. This resul t was 
also seen in semi-analytic work bv lPower et al.l ((20101) who 
found that cold gas was being turned into stars in these 
low-mass systems. 

We summarise the HI evolution in Fig. where we 
integrate the Schechter function described in Fig. [8] even 
though we can clearly see in those figures that the low-mass 
tail undergoes significant evolution in model ZC^WFB and 
is shallower than measured locally. The high-mass tail ap- 
pears to show much more gradual evolution over this period. 
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Figure 7. The HI, HII, H2 and stellar half-mass radii (in proper units) as a function of halo mass for 2 = 0, 1 and 2, left to right 
respectively. We consider runs with metal-line cooling and supernovae feedback, ZCJ\/VFB in the top row and also with the addition of 
AGN in ZC-WFB^AGN in the bottom row. Power law fits to all radii-mass relations are shown, with best-fit values given in Tabic [5] 



Observers will be constrained to observing only the largest 
systems so to, crudely, mime this selection we limit the fit 
of the Schechter function to HI masses above 10^" h~^MQ 
and fix the faint end to —1.5 (the value we find at z = 2). 
In this case we well agree with observations at low redshift 
and predict little evolution will be measured to high red- 
shift. We show a number of observational measurements in 
Fig. |9] and, while we did not attempt to make an exhaus- 
tive compilation, we consider this a reasonable summary of 
the current understanding of the cosmic HI density evo- 
lut ion. The results ra nge f rom direct, local H I detections 
bv lZwaan et all (|2005l ') and lMartin et all (|20ld ) to a stack- 
ing of the HI signal at a moderate redshift by iLah et al.l 
((20071). Then there are inferred cosmic HI densities through 
the use of HST observation s of low-redshift DLA systems 
traced by Mgll absorbers (jRao et al.l 120061: Meiring et al.l 
boill) and then at higher red s hift b y iProchaska fc Wolfg 
(|2009h and iNoterdaeme et all (|2009l ') using SDSS quasar 
spectra. There is disagreement between these last two re- 
sults due, primarily, to a more aggressive treatment for the 
bias affecting strong absorbers in the latter work. 



5.2 Molecular and Cold Gas Mass Functions 

We now consider the molecular mass function, shown in the 
left column of Fig. IIOI and the more general cold gas mass 
function, shown in the right column. Cold gas is defined as 
all gas that is gravitationally bound to the halo and that 
has temperature T < 10^'^K, as well as all star-forming gas. 
Note that neither the cold or molecular gas mass depends 
on the self-shielding measure used. 

For the z — molecular hydrogen mass function, top 
left panel, we have only a limited mass range (Mhj — 
1 — 2 X 10^° H'^Mq), but the abundance of objects is in 
reasonably good agreement with the data of iKeres et al.l 
((2003), using a constant (X-factor) conversion between CO 



and H2 mass. (We implicitly used t his conversion rather than 
the va riable measure advocate d bvlObreschkow fc Rawlingd 
l2009bl . when we adopted the iLerov et al.l 120081 molecular- 
atomic abundance normalisation.) It is clear that the high- 
mass end of the Schechter function evolves strongly as a 
function of redshift, with many more high molecular mass 
systems at higher redshift. The overall cosmic molecular 
mass density at 2: = 2 is predicted to be a factor two larger 
than observed locally. Thus, the large increase in stellar den- 
sity between z = 2 and the present day is reflected in the 
molecular mass function rather than the HI mass function, 
which appears to be demonstrate mild evolution. This over- 
all effect is also seen for a range of other physics schemes 
probed in Fig. 1121 

The cold gas mass function can be estimated observa- 
tionally by converting the HI mass function from ALFALFA, 
corrected for the He fraction assuming a molecular - atomic 
ratio based on the ratio o f the global densities of each. As ad- 
vocated by I Power et al.l (|20ld ). for a flxed molecular-atomic 
ratio this conversion is akin to dividing the HI mass by 0.54. 
For a conversion factor that depends on galactic scale p rop- 
erties, as considered in lObreschkow fc Rawlingd (|2009bl ) the 
HI mass - cold mass co nversion factor is given in Equation 
4. in lPower et"all (|2010l ') as 



MHi = 0.76Mcoid/(l + i?s. 



(10) 



where we previously defined -Rsurf in Equation [l] The exact 
conversion of HI mass to cold gas mass is especially impor- 
tant at the high-mass end. If we use a variable conversion 
ratio, then our simulated gas masses are too low at the high- 
mass end. The situation is reversed in the case of a constant 
conver sion factor. For the case of the constant conversion 
factor, I Power et al.l (|2010l ) find that the semi-analytic pre- 
scriptions also typically overproduce the number of gas-rich 
systems. There appears to be a generic disagreement be- 
tween both hydro and semi-analytic simulations with the 
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Figure 8. The HI mass function from ZC-WFB at 2 = 0, 1 and 2; top, middle and bottom respectively, with the 2 = ALFALFA HI 
mass function shown for comparison. Each column displays results for the different self-shielding methods: the left column is method 
SSNone (no self-shiclding); the middle column is method SSalfalfa (pressure threshold tuned to match the ALFALFA data in the 
ra nge Mhi = lO^'^ — 10^^ h~ ^M.p, ) and the right column is method SSdla (pressure threshold tuned to match the 2 = 2, Odla results 
of lProchaska fe Wolfj 120091) . Typically, the pressure value found to match the 2 = 2 DLA results is higher than that used for 2 = 0, 
allowing us to bracket the available range of possible HI mass functions for the ZC-WFB simulations. Poisson errors are assumed for 
individual simulation data points. We also show a Schechter function fit to the simulated data (for methods SSalfalfa a-nd SSdla) 
and present best-fit values for a, 8 and M* (and their errors, from bootstrap resampling) in the legend. Values for Hhi, calculated from 
the integrated Schechter function, are also shown. Clearly, rijji — S^qla for method SSqla at 2 = 2 but Qhi < ^ALFALFA in method 
SSalfalfa (i-c this model predicts less HI at high-redshift). 



data if a constant HI/H2 is assumed. This lends increasing 
support to the idea that the ratio depends on the global 
properties of the system in which it is meas ured, e.g. as con- 
sidered in lObreschkow fc Rawlingj (|2009a ). 



EFFECTS OF FEEDBACK AND METAL 
COOLING ON HI MASS FUNCTION 



As w e have demonstrated in previous work (JDuffv et al.l 
I2OI0I ). the gas distribution in a halo can vary significantly 
with the physics implementation. This is especially true as 
the ZC_ WFB model is unable to remove gas from the largest 
galaxies, and hence is unable to form red ellipticals by sup- 
pressing star formation in the most massive systems. We 
attempt to analyse the sensitivity of our results to the wide 
array of possible feedback types, as well as cooling schemes, 



used in this previous work. We have created the HI mass 
function at z = and 2 for the different physics prescrip- 
tions and present these in Fig. Illl again applying the HI 
mass cut from Equation [H] and tuning each physics simula- 
tion independently in the same manner as Section [3. II with 
pressure values given in Table O Equivalent plots to Fig. IIUI 
for the H2 and cold gas mass functions, are shown in Fig. 1121 

Without any self-shielding implemented at 2; = (top- 
left panel of Fig. Ilip the simulations with strong feedback 
{ZC.WFBuiGN and ZC.SFB) result in a greater num- 
ber of massive HI systems than the schemes without effec- 
tive feedback. This is a result of more cold gas being pre- 
vented from forming stars, as seen in the top right panel 
of Fig. 1121 When we tune the self-shielding prescription to 
match the ALFALFA HI cosmic density (over the HI mass 
range lO^'' — 10^^ /i~^Mq; method SSalfalfa) the difference 
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Figure 10. The simulated H2 (cold gas) mass function is shown in th e left (right) colum n at ^ = (top), 1 (middle) and 2 (bottom), 
with the z = H2 mass function r esults from the CO mass f unction of lKeres et alj 1120031 ) when using a constant (X-factor) conversion 
or the varying scheme proposed bv lObreschkow fc Rawlingd l|2009lJ ) (denoted with the dotted and dashed lines, respectively). The cold 
gas mass function is estimated by converting the ALFALFA detections through either a c onstant assumed at omic-molecular ratio or a 
variable scheme that depends on the disc gas and stellar mass, with both values taken from lPower et alJ l|2010l '). The errors on the points 
are Poissonian. We have made a resolution cut in Mh2 (A/cas) using Equation |6] Note that the molecular mass function evolves more 
strongly than the cold gas mass as the high mass end. 



between the various physics schemes at 2 = 0, top middle 
panel, is reduced as expected. 

At 2 = 2 we can tune to the DLA result (method 
SSdla), shown in the bottom right panel, which also shows 
very little difference between the schemes, although the run 
with no feedback (PrimC'-NFB) shows a lack of massive HI 
systems in comparison to the other runs. We see from the 
molecular and cold gas mass functions at 2 = 2 (bottom 
panels of Fig. I12|) that more molecular hydrogen is present 
at high-redshift, due to both the lack of feedback and time 
to use up the gas and form stars. 

It appears that the HI mass function is at least as sen- 
sitive, if not more so, to the modelling of self-shielding as 
to the sub-grid baryonic physics. This could be seen as an 
unexpected result and highlights the need for a more careful 
modelling of the optically thick gas in and around galaxies. 

There is significantly more variation in the molecular 



hydrogen mass function when the baryonic physics is varied, 
as shown in the left column of Fig. [12] The situation at z = 2, 
bottom left panel, is complex but apparently the strong feed- 
back schemes remove large amounts of the densest gas rel- 
ative to the simulations without feedback. The cold gas at 
Mcas JS 10®'^ h~^MQ is, on the other hand, completely in- 
sensitive to the sub-grid physics scheme used provided there 
is some feedback to suppress the rapid overcooling demon- 
strated by the PrimC^NFB model in the bottom right panel. 
The potential paradox presented here between the greater 
sensitivity of the H2 to the feedback schemes than the HI 
when the cold gas also appears insensitive is a consequence 
of the cold gas having no dependancy on the self-shielding 
prescription used to calculate the HI. As a result the HI 
fraction of the cold gas can vary to match the observations, 
giving a greater freedom in the HI results compared to the 
H2 which have no such variability. Thus the cold gas mass 
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Figure 11. As Fig. [S] but now showing a comparison of HI mass functions for runs with different physics prescriptions, at z = (top 
panels) and z = 2 (bottom panels). Again, the columns show results for self-shielding methods SSNonei SSalfAlfAi SSdlAi Isft to right 
respectively. The low-redshift result from ALFALFA is also shown for comparison. Note that certain simulations, e.g. PrimC^WFB and 
PrimC-NFB, have too little cold gas at low redshift that even when setting the Pshield = 0, converting all cold gas in the galaxy to be 
neutral. This is seen most clearly in the top-middle panel in which the red dashed and blue dotted lines lie below the other simulations 
at all mass bins. 
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Figure 9. A compilation of some recent measurements of the 
cosmic HI densities at various epochs. We have included the in- 
tegrated values for the simulation ZC.WFB (with self-shielding 
method SSdla) to demonstrate the difficulty in reproducing the 
observed HI density. Even though at z = we can reproduce 
the mass function about the 'knee' of the Schechter function, wc 
do not have enough low-mass objects to recover the observed 
Ohi- However, if we fix the faint end slope to the value at 2 = 2 
(a = -1.5) and limit the fit to haloes with Mhi > 10^° /i'^Mq 
then we recover the blue diamond points which are in better 
agreement with the low redshift results. They also demonstrate 
that, in the simulations at least, the evolution in HI is con- 
strained to th e faint end. Note that we h ave tuned the 2 = 2 case 
to agree with [ P rochas ka fc Wolfd 1 120091 ) rather than the results 
of iNoterdaeme et al.l 1120091 ). see text. In our models the evolu- 
tion of Ofil is driven by the low-mass tail, whereas most 21cm 
observations will be unable to detect this population. 



across simulations is similar but the distribution of this gas 
as a function of pressure is not, which is then most easily 
visible in the H2 mass functions shown in Fig. 1121 



7 CONCLUSIONS 

We have used hydrodynamical simulations from the OWLS 
project to investigate the distribution of neutral (atomic HI 
and molecular H2) hydrogen in cosmological volumes across 
the redshift range z — — 2. The large suite of simula- 
tions allowed us to probe both the effects of numerical res- 
olution and the sub-grid baryonic physics implemented (es- 
pecially the effects of metal enrichment on the gas cooling 
rate and the different methods and sources of feedback). In 
addition, we considered several methods to correct for the 
lack of self-shielding in our simulations in moderate-high- 
density regions, where the assumption that the gas is opti- 
cally thin to a uniform UV/X-ray background breaks down. 
With this in mind, we have found that, together with an 
empirical pressure-based law to discern the molecular frac- 
tion in galaxies, the distribution of HI can be successfully 
simulated in our cosmological volumes. We have succeeded 
in recreating the observed local HI mass function for ob- 
jects with HI mass greater than 10^'' /i^'^M© but with a 
significantly shallower faint-end slope. We consider evolu- 
tion across the redshift range z = — 2 and detect signifi- 
cant increases in low-mass HI objects which corresponds into 
a small, yet systematic, increase in the HI cosmic density 



with increasing redshift. If we mime an observational mag- 
nitude limited survey, namely objects with HI masses above 
10^" /i~^Mq, and fix the faint end slope to the 2 = 2 value 
then we find that the inferred cosmic HI density is in better 
agreement with observations at 2: = and that observers 
will infer little evolution in the period 2 = — 2. During this 
time the molecular hydrogen component evolves strongly, 
with many more massive systems in the distant past than 
now, in qualitative agreement with lObreschkow fc Rawlingd 
(|2009al ). 

When comparing the predictions of the simulations with 
varying physics prescriptions, including an AGN feedback 
model, we found that the HI mass functions were insensi- 
tive to these treatments. Provided one can tune the pressure 
based self-shielding threshold, we can recover the overall HI 
distribution to within a factor of two using just a single 
threshold value. 

However, the implementation of the self-shielding 
threshold has at least as large an effect on the HI mass 
function than the differences between the sub-grid cooling 
and feedback models we consider here, these models essen- 
tially bracket the extrema of sub-grid implementations and 
hence to improve on this level of accuracy in t he future we 
will u tilise radiative transfer schem es such as in lAltav et al.l 
(|2010l ) and lPawhk fc Schavel (|2008l ). 

For the case of the simulation with supernovae feedback 
and realistic cooling we found a peak HI - halo mass fraction 
at a given halo mass across our redshift range, 2 = — 
2, which we attribute to two effects. The first is that in 
more massive haloes the neutral component is dominated 
by denser H2 while at lower masses the overall gas fraction 
of the halo has been reduced due to feedback. This peak is 
also found to be the 'knee' of the HI mass function, which is 
well reproduced in our simulation with supernova feedback. 

The typical sizes, denoted by the half-mass radius, of 
the HI, H2 and stars in this work were found to vary as a 
function of total halo mass and redshift. We found, how- 
ever that the HII radius increased with the cube root of the 
halo mass, characteristic of a smoothly distributed isother- 
mal profile. The typical sizes of the HI disk increased with 
redshift but that the steepness of this relation with halo 
mass decreased from a linear dependency to approximately 
a cube root (tracing of the more diffuse, ionised component). 
The stars and molecular hydrogen closely traced each other 
throughout redshift, with an identical half mass radius at 
galaxy scales, i.e. the pivot mass of 2 x 10^^ h~^MQ, but 
differed significantly at 2 = for galaxy group and cluster- 
sizes, with the stellar radius a factor 5 more extended the 
molecular component. We also tested the feedback from an 
AGN and found that the molecular and stellar components 
were more extended than their counterparts in the simula- 
tion without an AGN, at all redshifts. This is in agreement 
with I Sales et al.l (|2010f ) who studied galaxies at 2 = 2 and 
found that such strong feedback resulted in the formation 
of large, extended galaxies. 

Finally, we have found that the ratio of HI (H2) to stel- 
lar mass of our objects is negatively (positively) correlated 
with the stellar mass. The normalisation of the model with 
SNe feedback was found to be in particularly good agree- 
ment with observations. Furthermore, we predict that the 
ratio of HI (H2) to stellar mass at a given stellar mass 
increases significantly with redshift. This means that the 
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hydrogen signal from small systems will be easier to find 
than an extrapolation of from larger optical counterparts 
would suggest. This improves the chances that optically faint 
galaxies may be detectable in upcoming ultra-deep 21 cm 
and molecular hydrogen observations. 



ACKNOWLEDGEMENTS 

ARD was supported under an RAS Centre of Excellence 
grant and further acknowledges this work was begun during 
an STFC PhD studentship. The authors wish to thank Mar- 
tin Meyer, Rob Grain and Max Pettini for helpful and stim- 
ulating discussions. Further thanks to Kirsten Gottschalk, 
Serena Bertone and Freeke van de Voort for extensive tech- 
nical help. Additionally ARD would like to thank Attila 
Popping and Ann Martin for making available comparison 
data and useful discussions of results. The simulations pre- 
sented here were run on Stella, the LOFAR BlueGene/L 
system in Groningen and on the Cosmology Machine at 
the Institute for Computational Cosmology in Durham as 
part of the Virgo Consortium research programme. This 
work was sponsored by the National Computing Facilities 
Foundation (NCF) for the use of supercomputer facilities, 
with financial support from the Netherlands Organization 
for Scientific Research (NWO). This work was supported by 
a Marie Curie Initial training Network CosmoComp (PITN- 
GA-2009-238536). 



REFERENCES 

Abdalla F. B., Rawhngs S., 2005, MNRAS, 360, 27 

Altay G., Theuns T., Schaye J., Crighton N. H. M., Dalla 
Vecchia C, 2011, ApJ, 737, L37 

Bandara K., Crampton D., Simard L., 2009, ApJ, 704, 1135 

Barnes D. G., Staveley-Smith L., de Blok W. J. G., Oost- 
erloo T., Stewart I. M., et al., 2001, MNRAS, 322, 486 

Bate M. R., Burkert A., 1997, MNRAS, 288, 1060 

Booth C. M., Schaye J., 2009, MNRAS, 398, 53 

Booth C. M., Schaye J., 2010, MNRAS, 405, LI 

Chabrier C, 2003, PASP, 115, 763 

Chang T., Pen U., Bandura K., Peterson J. B., 2010, Na- 
ture, 466, 463 

Cole S., Lacey C, 1996, MNRAS, 281, 716 

Grain R. A., Eke V. R., Frenk C. S., et al., 2007, MNRAS, 
377, 41 

Dalla Vecchia C, Schaye J., 2008, MNRAS, 387, 1431 

Davis M., Efstathiou C, Frenk C. S., White S. D. M., 1985, 
ApJ, 292, 371 

Dolag K., Borgani S., Murante C, Springel V., 2009, MN- 
RAS, 399, 497 

Drory N., Bender R., Feulner G., Hopp U., Maraston C., 
Snigula J., Hill G. J., 2003, ApJ, 595, 698 

Drory N., Salvato M., Gabasch A., Bender R., Hopp U., 
Feulner G., Pannella M., 2005, ApJL, 619, L131 

Duffy A. R., Battye R. A., Davies R. D., Moss A., Wilkin- 
son P. N., 2008, MNRAS, 383, 150 

Duffy A. R., Moss A., Staveley-Smith L., 2011, ArXiv e- 
prints 

Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C, Battye 
R. A., Booth C. M., 2010, MNRAS, 405, 2161 



Erb D. K., Shapley A. E., Pettini M., Steidel C. C, Reddy 

N. A., Adelberger K. L., 2006, ApJ, 644, 813 
Fabello S., Catinella B., Giovanelli R., Kauffmann G., 

Haynes M. P., Heckman T. M., Schiminovich D., 2011, 

MNRAS, 411, 993 
Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., 

Kingdon J. B., Verner E. M., 1998, PASP, 110, 761 
Fontana A., Pozzetti L., Donnarumma I., Renzini A., 

Cimatti A., et al., 2004, A&A, 424, 23 
Giovanelh R., Haynes M. P., Kent B. R., et al., 2005, AJ, 

130, 2598 
Haardt F., Madau P., 2001, in Neumann D. M., Tran 

J. T. v., eds, Clusters of Galaxies and the High Red- 
shift Universe Observed in X-rays Modelling the UV/X- 

ray cosmic background with CUBA 
Haring N., Rix H., 2004, ApJL, 604, L89 
Kennicutt Jr. R. C, 1998, ARA&A, 36, 189 
Keres D., Yun M. S., Young J. S., 2003, ApJ, 582, 659 
Komatsu E., Smith K. M., Dunkley J., et al., 2011, ApJS, 

192, 18 
Krumholz M. R., McKee C. F., Tumhnson J., 2008, ApJ, 

689, 865 
Krumholz M. R., McKee C. F., Tumhnson J., 2009, ApJ, 

693, 216 
Lacey C, Cole S., 1994, MNRAS, 271, 676 
Lagos C. d. P., Baugh C. M., Lacey C. C, Benson A. J., 

Kim H.-S., Power C., 2011, ArXiv e-prints 
Lah P., Chengalur J. N., Briggs F. H., et al., 2007, MNRAS, 

376, 1357 
Lang R. H., Boyce P. J., Kilborn V. A., et al., 2003, MN- 
RAS, 342, 738 
Leroy A. K., Walter F., Brinks E., et al., 2008, AJ, 136, 

2782 
Madau P., Ferguson H. C, Dickinson M. E., Giavalisco M., 

Steidel C. C, Fruchter A., 1996, MNRAS, 283, 1388 
Martin A. M., Papastergis E., Giovanelli R., Haynes M. P., 

Springob C. M., Stierwalt S., 2010, ApJ, 723, 1359-1374 
McCarthy I. G., Schaye J., Ponman T. J., Bower R. G., 

Booth C. M., Vecchia C. D., Grain R. A., Springel V., 

Theuns T., Wiersma R. P. C, 2010, MNRAS, pp 740— f 
McQuinn M., Oh S. P., Faucher-Giguere C.-A., 2011, ArXiv 

e-prints 
Meurer G. R., Hanish D. J., Ferguson H. C, Knezek P. M., 

Kilborn V. A., et al., 2006, ApJS, 165, 307 
Meyer M. J., Zwaan M. A., Webster R. L., Staveley-Smith 

L., Ryan- Weber E., et al., 2004, MNRAS, 350, 1195 
Meiring, J. D., Tripp, T. M., Prochaska, J. X., Tumlinson, 

J., Werk, J., Jenkins, E. B., Thom, C, O'Meara, J. M., 

Sembach, K. R., 2011, ApJ, 732, 35 
Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319 
Noterdaeme P., Petitjean P., Ledoux C., Srianand R., 2009, 

A&A, 505, 1087 
Obreschkow D., Croton D., DeLucia G., Khochfar S., Rawl- 

ings S., 2009, ApJ, 698, 1467 
Obreschkow D., Rawlings S., 2009a, ApJL, 696, L129 
Obreschkow D., Rawlings S., 2009b, MNRAS, 394, 1857 
Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651 
Pawlik A. H., Schaye J., 2009, MNRAS, 396, L46 
Pontzen A., Governato F., Pettini M., et al., 2008, MN- 
RAS, 390, 1349 
Popping A., Dave R., Braun R., Oppenheimer B. D., 2009, 

A&A, 504, 15 



20 A. R. Duffy et al. 



Power C, Baugh C. M., Lacey C. G., 2010, MNRAS, 406, 

43 
Pozzetti L., Cimatti A., Zamorani G., Daddi E., Menci N., 

et al, 2003, A&A, 402, 837 
Prochaska J. X., Herbert-Fort S., Wolfe A. M., 2005, ApJ, 

635, 123 
Prochaska J. X., O'Meara J. M., Worseck G., 2010, ApJ, 

718, 392 
Prochaska J. X., Wolfe A. M., 2009, ApJ, 696, 1543 
Quinn T., Katz N., Efstathiou G., 1996, MNRAS, 278, L49 
Rao S. M., Turnshek D. A., Nestor D. B., 2006, ApJ, 636, 

610 
Razoumov A. O., Norman M. L., Prochaska J. X., Sommer- 

Larsen J., Wolfe A. M., Yang Y.-J., 2008, ApJ, 683, 149 
Sales L. V., Navarro J. F., Schaye J., Dalla Vecchia C, 

Springe! V., Booth C. M., 2010, MNRAS, 409, 1541 
Savagho S., Glazebrook K., Le Borgne D., Juneau S., Abra- 
ham R. G., et al., 2005, ApJ, 635, 260 
Schaye J., 2001a, ApJL, 562, L95 
Schaye J., 2001b, ApJ, 559, 507 
Schaye J., 2004, ApJ, 609, 667 

Schaye J., DaUa Vecchia C., 2008, MNRAS, 383, 1210 
Schaye J., Dalla Vecchia C., Booth C. M., Wiersma 

R. P. C., Theuns T., Haas M. R., Bertone S., Duffy A. R., 

McCarthy I. G., van de Voort F., 2010, MNRAS, 402, 

1536 
Seljak U., Zaldarriaga M., 1996, ApJ, 469, 437 
Spergel D. N., Bean R., Dore O., et al., 2007, ApJS, 170, 

377 
Springe! V., 2005, MNRAS, 364, 1105 
Tacconi L. J., Genze! R., Neri R., Cox P., Cooper M. C., 

et a!., 2010, Nature, 463, 781 
Tescari E., Vie! M., Tornatore L., Borgani S., 2009, MN- 
RAS, 397, 411 
Tremaine S., Gebhardt K., Bender R., Bower G., Dressier 

A., et a!., 2002, ApJ, 574, 740 
Tremonti C. A., Heckman T. M., Kauffmann G., Brinch- 

mann J., Chariot S., et al., 2004, ApJ, 613, 898 
van de Voort F., Sdiaye J., Booth C. M., Dalla Vecchia C, 

2011, MNRAS, pp 805-f 
van de Voort F., Schaye J., Booth C. M., Haas M. R., DaUa 

Vecchia C, 2011, MNRAS, 414, 2458 
Wiersma R. P. C, Schaye J., Smith B. D., 2009, MNRAS, 

393, 99 
Wiersma R. P. C, Schaye J., Theuns T., DaUa Vecchia C, 

Tornatore L., 2009, MNRAS, 399, 574 
Zwaan M. A., Meyer M. J., Staveley-Smith L., Webster 

R. L., 2005, MNRAS, 359, L30 
Zwaan M. A., Staveley-Smith L., Koribalski B. S., et al., 

2003, AJ, 125, 2842 



APPENDIX A: RESOLUTION TESTS 
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Figure Al. Resolution tests for the HI mass in the FoF group as 
a function of M'?^'" in ZC- WFB. The top panel shows results for 
the 50 h^^ Mpc box at z = 0, while the bottom panel is for the 
25 h^'^Mpc box at z = 2. At z = 2 (0) we show results for 3 (2) 
simulations of the same box, spanning a factor of 64 (8) in particle 
number/mass. The vertical lines indicate the mass limit given 
in Equation [6] (roughly speaking this is equivalent to 7,000 dark 
matter particles) for each run. The curves correspond to haloes an 
order of magnitude in mass below this limit; the results appear 
converged, therefore we can be confident with our conservative 
limit in Equation [6] We bin the haloes that have masses greater 
than the vertical mass limits and denote their median values with 
symbols (triangles for the solid curve, squares for the dotted). 
Errors are la values, estimated by bootstrap analysis and added 
in quadrature with Poissonian error estimates. 



the particular baryonic species under consideration) and the 
particle mass. We also show that our results cannot be sig- 
nificantly affected by the lack of structure on scales larger 
than our available box sizes. Although we consider the case 
of the neutral phases of hydrogen, similar tests have been 
performed for the stars and ionised hydrogen to create the 
mass limits quoted in Equation |6] 



Here we investigate how numerical resolution and simula- 
tion box size affect our main results. In particular, we study 
how increasing the number of particles (and decreasing the 
gravitational softening length) affects the properties of the 
same sample of haloes. We also investigate the effects of 
increasing the simulation volume while keeping the reso- 
lution constant. As a result, we identify a simple scaling 
between the minimum required halo mass (that depends on 



Al HI mass range 

In order to determine the minimum halo size, and the corre- 
sponding HI mass resolution, we select FoF groups from the 
' ZC-WFB' simulation according to their total mass, M^j''". 
After grouping the objects in fixed mass bins we calculate 
the median HI mass of the bin. The variation in median 
HI mass with halo mass is shown for two simulation sets in 
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Figure A2. As Fig. lAll but comparing different box sizes (25, 
50 & 100 /i~^Mpc, where available) at fixed resolution (where 
available) at ^ = and 2. 



Fig. lAll with increasing resolution in the 50 h~^Mpc box at 
z = (top panel); and those with increasing resolution in 
the 25/i~^Mpc box at z = 2 (bottom panel). 

As one would expect, the curves demonstrate a positive 
correlation between HI and halo mass. Those haloes which 
lie on the same relation between simulations typically have 
at least ~ 7, 000 dark matter particles. We adopt a conserva- 
tive mass limit, shown as a vertical line for each simulation, 
and identify the HI mass corresponding to this halo mass 
for each simulation. As mentioned in Section [4. II we found 
that we could summarise this limit for all species as a simple 
power-law scaling relation of the simulation mass resolution, 
as given in Equation. [6] in the main body of the text. In our 
study, we have applied this limit to the baryonic masses, 
resulting in an even more conservative lower HI mass limit. 

We also compare the same relation for different sim- 
ulation volumes at the same resolution in Fig. IA2I These 
demonstrate that the results are insensitive to the total vol- 
ume simulated. 



A2 HI mass function 

Given our established lower HI mass limits, we now investi- 
gate the range of masses over which we can trust the HI 
mass function. To that end we show the HI mass func- 
tions for the simulations (going down an order of magni- 
tude in mass) in Fig. IA3I illustrating our minimum HI mass 
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Figure A3. Effects of resolution on the HI mass function for the 
50 fe^^Mpc boxes at z = and the 25 h~^Mpc boxes at ^ = 2 
(top and bottom rows respectively). These results span a factor 
of 64 in particle mass between the lowest and highest resolutions. 
Lower HI mass limits are shown as vertical linos, above which 
the abundance of haloes is reasonably well converged (Poisson 
errors in abundance are assumed). The mass function starts to 
rises more steeply below the resolution limit, suggesting a nu- 
mcri cal artefact. We have also included the z = data-points 
from I Popping et al.l 1120091 ) in the top panel. We note that our HI 
mass limit applied to their simulation (orange line) suggests that 
their results may suffer from numerical effects. 



cut-off as a vertical line. For comparison we also overplot 
the datapoints from P09; the resolution of their simulation 
(with 256"^ dark matter particles in a 32 h~^Mpc box) lies in 
between our 100/i^^Mpc and 50/i^^Mpc simulations with 
512"' dark matter particles. 

It is clear that our mass function is reasonably well con- 
verged above our HI mass limits; below these limits there is 
an artificial rise in the number of haloes, an effect that may 
be present in the results of P09 (also shown in the figure) 
who study the mass function down to 3.5 x 10* /i~^Mq; we 
would estimate their resolution limit to be cut at masses 40% 
higher, close to where their simulations exhibit the charac- 
teristic rapid rise in the HI mass function. We have also 
checked the effects of increasing volume at fixed resolution 
and find excellent convergence between the simulations. 



22 A. R. Duffy et al. 



A3 Molecular Hydrogen 

The molecular hydrogen in the simulations comes from 
higher density regions in galaxies than the HI and there- 
fore may be expected to be a more demanding quantity to 
resolve. Again, by comparing the molecular mass from sim- 
ulations with different resolution but of the same volume, 
we can determine a resolution-dependent mass limit. As ex- 
pected, our H2 mass limit is somewhat higher (by around 
an order of magnitude) than the HI mass limit, as given in 
Equation [6] 

In Fig. IA4I we show the molecular H2 mass functions 
for simulations with increasing resolution, in the 50 /i~^Mpc 
box at z — (top panel) and the 25/i~^Mpc box at 2: = 2 
(bottom panel). Vertical lines denote the resolution limits 
for the simulations. The overlapping curves demonstrate a 
reasonably well converged result, although it is clear that the 
mass resolution is extremely demanding. The z = results 
from P09 are also shown (we estimate their limit is 5 x 
10^ /i~^M0, they choose a limit 75 times less than this). 

Finally, we have also tested the sensitivity of the molec- 
ular mass function to the effects of large-scale structure by 
reducing the simulation volume at constant mass resolution. 
As with the HI case, there appears to be excellent agreement 
between the overlapping box sizes. 
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Figure A4. Molecular H2 mass functions for simulations witii 
different resolution. Tlie top panel is for the 50/i~^Mpc simu- 
lation box at z = and the bottom panel is for the smaller 
25/i~^Mpc box at z = 2. Resolution limits are shown as vertical 
lines and we assume a Poisson error for eac h data-point. We have 
also included the 2 = data-points from I Popping et al.l 1120091 ) 
in the top panel. As in Fig. IA3I our resolution limit applied to 
their simulation (orange line) suggests that their results may suf- 
fer from numerical effects. The grey curves are th e z = H2 
mass function results from the CO mass function of lKeres et al.l 
(I2OO3) when using a c onstant (X-factor) conver s ion or t he varying 
scheme proposed by lObreschkow fc Rawlingj l|2009b|) (denoted 
with the dotted and dashed lines, respectively). 



